!     -------------
!     PROGRAM PROG1
!     -------------
!
!     PLANETARY TABLES FROM 1900 TO 2100.
!
!     Ref : JC-GF February 2000.
!
! --- Object -----------------------------------------------------------
!
!     Example for using the subroutine PLATAB:
!     - Reading the ASCII planetary files.
!     - Computation of rectangular heliocentric coordinates of the
!       planets (equator and equinox of J2000).
!
!     First Julian date: JD2415020.5.
!     Step length: 14500.25 days.
!     Number of dates: 5.
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
      character*10 pla(9)
      character*12 fich
!
      dimension r(6)
!
      data pla/
     .'Mercury','Venus','E-M Baryc','Mars','Jupiter',
     .'Saturn','Uranus','Neptune','Pluto'/
!
! --- Parameters -------------------------------------------------------
!
      t0=2415020.5d0
      step=14500.25d0
      ndat=5
!     lu=10
!
! --- Results file (unit 20) -------------------------------------------
!
      fich='prog1.txt'
      open (20,file=fich)
      write (20,1001)
      write (*,1001)
!
! --- Compute the planetary coordinates --------------------------------
!
      do ipla=1,9
         t=t0
         do n=1,ndat
            call PLATAB (t,ipla,lu,r,ierr)
            if (ierr.ne.0) then
               write (*,1002) ierr
               stop
            endif
            write (20,1003) pla(ipla),t,r
            write (*,1003)  pla(ipla),t,r
            t=t+step
         enddo
      enddo
      write (20,1004) fich
      write (*,1004)  fich
      stop
!
! --- Formats ----------------------------------------------------------
!
1001  format (/2x,'PLANETARY TABLES : PROG1'/
     .        /2x,'Rectangular coordinates (equator and equinox J2000)')
1002  format (/2x,'ERREUR : ',i2/)
1003  format (/2x,a/2x,'JD:',f16.5/
     .       2x,'x :',f16.12,2x,'y :',f16.12,2x,'z :',f16.12,' au'/
     .       2x,'x'':',f16.12,2x,'y'':',f16.12,2x,'z'':',f16.12,' au/d')
1004  format (/2x,'DONE - Results File : ',a/)
!
      end
!
!
!
      subroutine PLATAB (tjd,ipla,lu,r,ierr)
! ----------------------------------------------------------------------
!
!     PLANETARY TABLES FROM 1900 TO 2100
!
!     Ref : JC-GF February 2000
!
! --- Object -----------------------------------------------------------
!
!     Compute planetary coordinates with tables built by a method of
!     representation using frequency analysis. For each planet there
!     is an ASCII file containing the corresponding planetary tables.
!
!     Coordinates: rectangular heliocentric coordinates.
!     Source: DE403 (Jet Propulsion Laboratory).
!     Frame: Dynamical Mean Equinox and Equator of J2000 (DE403).
!     Time: Barycentric Dynamical Time (TDB).
!
! --- Input ------------------------------------------------------------
!
!     tjd :       Julian date TDB (double real).
!
!     ipla :      Planet index (integer).
!                 1 : Mercury.
!                 2 : Venus.
!                 3 : Earth-Moon Barycenter.
!                 4 : Mars.
!                 5 : Jupiter.
!                 6 : Saturn.
!                 7 : Uranus.
!                 8 : Neptune.
!                 9 : Pluto.
!
!     lu :        Logical unit of the ASCII file (integer).
!
! --- Output -----------------------------------------------------------
!
!     r(6) :      Table of rectangular coordinates (double real).
!                 r(1) : X  position (au).
!                 r(2) : Y  position (au).
!                 r(3) : Z  position (au).
!                 r(4) : X' velocity (au/day).
!                 r(5) : Y' velovity (au/day).
!                 r(6) : Z' velocity (au/day).
!
!     ierr :      Error index (integer).
!                 0 : no error.
!                 1 : planet index error.
!                 2 : file error.
!                 3 : date error.
!
! --- Files ------------------------------------------------------------
!
!     tab1.mer : Tables for Mercury
!     tab2.ven : Tables for Venus
!     tab3.emb : Tables for Earth-Moon Barycenter
!     tab4.mar : Tables for Mars
!     tab5.jup : Tables for Jupiter
!     tab6.sat : Tables for Saturn
!     tab7.ura : Tables for Uranus
!     tab8.nep : Tables for Neptune
!     tab9.plu : Tables for Pluton
!
!     The files have to be present in the current directory.
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
      logical fexist
!
      character*1 num(9)
      character*3 nom(9)
      character*8 fname
!
      parameter (maxor=2,maxfq=250)
!
      dimension r(6)
      dimension v(3),w(3),ws(3)
      dimension nf(0:maxor),fq(maxfq,0:maxor)
      dimension sec(0:3,3)
      dimension ct(maxfq,0:maxor,3)
      dimension st(maxfq,0:maxor,3)
!
      data num/'1','2','3','4','5','6','7','8','9'/
      data nom/'mer','ven','emb','mar','jup','sat','ura','nep','plu'/
!
! --- Clear the results table ------------------------------------------
!
      do i=1,6
         r(i)=0.d0
      enddo
!
! --- Check up the planet index ----------------------------------------
!
      ierr=1
      if (ipla.lt.1.or.ipla.gt.9) return
!
! --- Open the ASCII file ----------------------------------------------
!
      fname='tab'//num(ipla)//'.'//nom(ipla)
      inquire (file=fname,exist=fexist)
!
      ierr=2
      if (.not.fexist) stop
!
      open (lu,file=fname,access='sequential',status='old',iostat=nerr)
      if (nerr.ne.0) return
!
      rewind lu
      read (lu,fmt='(12x,f10.1,f8.0,3i3)') tzero,dt,mx,imax,iblock
!
! --- Check up the date ------------------------------------------------
!
      ierr=3
      tmax=tzero+dt*iblock
      if (tjd.lt.tzero-0.5d0.or.tjd.gt.tmax+0.5d0) return
!
      nb=(tjd-tzero)/dt+1
      if (tjd.le.tzero) nb=1
      if (tjd.ge.tmax) nb=iblock
!
! --- Read the frequencies ---------------------------------------------
!
      do m=0,mx
         read (lu,fmt='(i4)') nf(m)
         do i=1,nf(m)
            read (lu,fmt='(f20.16)') fq(i,m)
         enddo
      enddo
!
! --- Read the separator -----------------------------------------------
!
      read (lu,*)
!
! --- Read the series --------------------------------------------------
!
      do k=1,nb
      do iv=1,3
!
!        Sub-interval and variable numbers
!
         read (lu,*)
!
!        Read the secular terms
!
         do i=0,imax,2
            read (lu,fmt='(2f14.0)') sec(i,iv),sec(i+1,iv)
         enddo
!
!        Read the Poisson terms
!
         do m=0,mx
            ip=mod(m,2)
            do i=1,nf(m)
               read (lu,fmt='(2f14.0)') a,b
               if (ip.eq.0) then
                  ct(i,m,iv)=a
                  st(i,m,iv)=b
               else
                  ct(i,m,iv)=b
                  st(i,m,iv)=a
               endif
            enddo
         enddo
!
      enddo
      enddo
!
! --- Time variables ---------------------------------------------------
!
      tdeb=tzero+(nb-1)*dt
      x=2.d0*(tjd-tdeb)/dt-1.d0
      fx=x*dt/2.d0
!
! --- Compute the positions (secular terms) ----------------------------
!
      do iv=1,3
         v(iv)=0.d0
         wx=1.d0
         max=2*imax-1
         do i=0,max
            v(iv)=v(iv)+sec(i,iv)*wx
            wx=wx*x
         enddo
      enddo
!
! --- Compute the positions (Poisson terms) ----------------------------
!
      wx=1.d0
      do m=0,mx
         nw=nf(m)
         do iv=1,3
            ws(iv)=0.d0
         enddo
         do i=1,nw
            f=fq(i,m)*fx
            cf=cos(f)
            sf=sin(f)
            do iv=1,3
               ws(iv)=ws(iv)+ct(i,m,iv)*cf+st(i,m,iv)*sf
            enddo
         enddo
         do iv=1,3
            v(iv)=v(iv)+ws(iv)*wx
         enddo
         wx=wx*x
      enddo
!
! --- Compute the velocities (secular terms) ---------------------------
!
      wt=2.d0/dt
      do iv=1,3
         w(iv)=0.d0
         wx=1.d0
         max=2*imax-1
         do i=1,max
            w(iv)=w(iv)+i*sec(i,iv)*wx
            wx=wx*x
         enddo
         w(iv)=wt*w(iv)
      enddo
!
! --- Compute the velocities (Poisson terms) ---------------------------
!
      wx=1.d0
      do m=0,mx
         nw=nf(m)
         do i=1,nw
            fw=fq(i,m)
            f=fw*fx
            cf=cos(f)
            sf=sin(f)
            do iv=1,3
               stw=st(i,m,iv)
               ctw=ct(i,m,iv)
               w(iv)=w(iv)+fw*(stw*cf-ctw*sf)*wx
               if (m.gt.0) w(iv)=w(iv)+m*wt*(ctw*cf+stw*sf)*wy
            enddo
         enddo
         wy=wx
         wx=wx*x
      enddo
!
! --- Stock the results ------------------------------------------------
!
      do iv=1,3
         r(iv)=v(iv)/1.d10
         r(iv+3)=w(iv)/1.d10
      enddo
!
! --- Exit -------------------------------------------------------------
!
      ierr=0
      close (lu)
      return
      end
